Algoritam je zasnovan na maloj Fermaovoj teoremi $x^{p - 1} \equiv 1 (mod$ $p)$, gde je $p$ prost broj
i pronalaženju netrivijalnih korena jedinice $x^2 \equiv 1 (mod$ $p)$ akko $p|(x^2 - 1) = (x - 1)(x + 1)$, tj. $x \equiv -1 (mod$ $p)$ ili $x \equiv 1 (mod$ $p)$.
Ako algoritam vrati vrednost $False$, broj je sigurno složen,
u suprotnom, ako je vraćena vrednost $True$, broj je verovatno prost sa verovatnoćom $p = 4^{-k}$, gde je $k$ parametar algoritma.
Detaljnije ovde
import random
# Pomocna funkcija za stepenovanje po modulu
def mod_pow(a, n, m):
result = 1
a = a % m
while n > 0:
if n % 2 == 1:
result = (result * a) % m
a = (a * a) % m
n = n // 2
return result
def miller_rabin(n, k):
if n <= 3:
if n == 1:
return False
return True
# n prost => n neparan => n = (2 ^ r) * d + 1
d = n - 1
r = 0
while d % 2 == 0:
r = r + 1
d = d // 2
for i in range(k):
a = random.randrange(2, n - 1)
x = mod_pow(a, d, n)
if x == 1 or x == n - 1:
continue
wittness = True
for j in range(r - 1):
x = mod_pow(x, 2, n)
if x == 1:
return False
if x == n - 1: # n - 1 = -1 (mod n)
wittness = False
break
if wittness:
return False
return True
k = 20
print(f'3 prost? {miller_rabin(3, k)}')
print(f'4 prost? {miller_rabin(4, k)}')
print(f'561 prost? {miller_rabin(561, k)}')
print(f'1021 prost? {miller_rabin(1021, k)}')